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ABSTRACT 


This  thesis  examines  two  implementations  of  the  para¬ 
bolic  equation  approximation  to  the  acoustic  wave  equation 
aimed  at  removing  three  errors  inherent  to  the  wide-angle 
parabolic  equation  (WAPE)  model.  First,  the  selection  of 
the  range-step  size  used  by  the  split-step  Fourier  algorithm 
affects  the  convergence  of  the  solution.  Second,  in  certain 
ocean  environments  WAPE  incorrectly  computes  the  down-range 
transmission  loss.  Finally,  WAPE  does  not  reproduce  the 
standard  normal  mode  basis  set  as  defined  by  normal  mode  the¬ 
ory.  A  double-precision  implementation  of  the  WAPE  (DP- 
WAPE)  is  developed  to  evaluate  the  dependence  of  solution 
convergence  on  the  numerical  precision  of  the  model . 

Finally,  an  implementation  that  is  insensitive  to  the  choice 
of  the  reference  sound  speed  (COIPE)  is  evaluated  for  its 
ability  to  reduce  or  remove  the  latter  two  of  these  three 
errors.  The  stability  of  the  WAPE  solution  was  found  to  be 
unaffected  by  the  DP-WAPE  implementation.  The  range-step 
dependence  is  inherent  to  the  split-step  algorithm.  The 
COIPE  corrects  the  transmission  loss  anomaly  and  satisfacto¬ 
rily  reproduces  the  standard  normal  mode  basis  set. 
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I.  INTRODUCTION 


The  fall  of  the  former  Soviet  Union  and  the  resulting 
shift  in  the  international  political  power  base  has  resulted 
in  a  dramatic  shift  in  both  the  content  and  the  application 
of  United  States  Naval  forces  during  the  past  decade.  The 
world's  largest  blue  water  navy  finds  itself  increasingly 
involved  in  green  water  operations.  This  shift  to  littoral 
waters  has  brought  to  the  forefront  the  ever  present  danger 
of  easily  obtainable  and  extremely  capable  diesel  submarines 
and  shallow  water  mines.  Countries  with  otherwise  limited 
military  resources  can  acquire  the  ability  to  thwart  or  at 
the  least  hinder  the  application  of  national  objectives.  As 
a  result  of  these  changes,  numerous  efforts  toward  the 
development  of  more  accurate  acoustic  propagation  prediction 
models  have  arisen.  These  include  areas  such  as  matched- 
field  processing,  transient  localization,  and  underwater 
acoustic  communications.  Essential  to  these  efforts  is  the 
development  of  acoustic  prediction  models  which  can  produce 
both  valid  and  accurate  results  in  highly  variable  shallow 
water  ocean  environments . 
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A.  OCEAN  ACOUSTIC  MODELING  TECHNIQUES 

Current  acoustic  modeling  efforts  generally  fall  into 
one  of  four  broad  categories.  These  are  ray  models,  wave- 
number  integration  techniques,  normal  mode  models,  and  para¬ 
bolic  equation  models.  Each  of  these  are  based  on  a  variety 
of  approximations  to  the  linear  acoustic  wave  equation. 

Many  of  the  early  models  are  based  on  the  geometrical  limit 
and  define  ray  trajectories  of  sound  propagation.  The 
resulting  ray  models  can  quickly  predict  pulse  propagation 
travel  times  and  provide  rough  approximations  of  the  sound 
field.  The  limitation  of  this  approach  is  its  reliance  on  a 
high-frequency  limit  which  neglects  finite- frequency  effects 
such  as  diffraction.  In  addition,  it  has  been  shown  (Smith, 
et.  al.,  1992)  that  the  ray  equations  form  a  set  of  coupled, 
nonlinear  equations  which  suffer  from  chaotic  solutions  in 
the  presence  of  range - dependenc e .  Given  the  current  impor¬ 

tance  of  low-frequency  propagation  in  range -dependent  envi¬ 
ronments,  the  usefulness  of  ray  models  in  shallow  water 
media  is  limited. 

Wavenumber  integration  methods  produce  highly  accurate, 
full-wave  solutions  to  the  wave  equation.  However,  they  are 
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computationaly  intensive  and  are,  by  design,  limited  to 
range- independent  environments.  They  can  be  quite  useful  in 
determining  near-field  effects  but  are  not  of  practical 
utility  in  predictions  of  far- field  acoustic  propagation. 

Normal  mode  models  are  based  on  a  separation  of  vari¬ 
ables  in  range  and  depth  of  the  acoustic  wave  equation. 

This  approach  maintains  the  full-wave  characteristics  of  the 
acoustic  field,  and  provides  highly  accurate  results  to  both 
simple  and  complex  range  dependent  ocean  environments .  The 
solutions  produced  are  based  on  the  single- frequency  Helm¬ 
holtz  wave  equation.  Due  to  the  computational  complexity 
involved  with  high  frequencies  and  broadband  pulses,  as  well 
as  range-dependent  media,  such  models  are  best  suited  for 
low-frequency,  shallow  water  environments  which  support  only 
a  small  number  of  propagating  modes.  The  results  remain  very 
accurate,  but  computational  time  can  become  counter  produc¬ 
tive  . 

The  final  acoustic  modeling  technique,  and  undoubtedly 
the  most  popular  in  range -dependent  environments,  is  based 
on  the  parabolic  approximation  to  the  Helmholtz  wave  equa¬ 
tion.  Numerous  marching  algorithms  have  been  developed  that 
solve  the  acoustic  field  as  a  boundary  value  problem.  While 
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it  is  a  one-way  wave  equation  treatment,  it  is  a  full-wave 
model  and  therefore  includes  all  of  the  finite- frequency 
effects  of  diffraction.  Two-way  adaptations  of  parabolic- 
equation  models  have  been  developed,  but  do  not  result  in 
significant  changes  for  the  majority  of  shallow  water  envi¬ 
ronments  of  interest.  Solutions  are,  however,  still  single 
frequency,  requiring  multiple  solution  generation  for  pulse 
propagation  responses.  It  is  this  method  that  will  be  the 
primary  focus  of  this  thesis. 

B.  HISTORY  OF  THE  PARABOLIC  EQUATION  APPROXIMATION 

The  parabolic  equation  (PE)  approximation  method  was 
first  applied  to  the  problem  of  underwater  acoustics  by  Tap- 
pert  (1977) .  Its  historical  roots  are  much  deeper,  however. 
Breakthroughs  in  physics  and  mathematics  in  the  mid-1920's 
provided  the  basis  for  the  parabolic  wave  equation  used  in 
acoustics  today.  The  standard  parabolic  equation  (SPE)  has 
the  same  form  as  the  Schrodinger ' s  equation  in  quantum 
mechanics.  The  earliest  documented  use  of  the  SPE  as  an 
approximation  to  the  theory  of  wave  propagation  dates  to  the 
mid-1940's  and  the  work  of  Leontovich  and  Fock  (1946)  who 
originally  coined  the  term  "parabolic  equation  method" . 


4 


They  applied  the  method  to  predicting  the  diffraction  on 
long  range  tropospheric  wave  propagation  caused  by  the 
spherical  shape  of  the  earth.  Fock  (1965)  expanded  this 
approach  to  include  high-frequency  scattering  as  well  as 
microwave  propagation  in  waveguides. 

With  the  advent  of  lasers  and  their  coherent  radiation 
sources  in  the  1960's,  the  PE  method  was  further  extended  to 
laser  beam  propagation  where  it  is  generally  referred  to  as 
the  "quasi-optical"  equation.  The  quasi-optical  equation  is 
most  widely  used  in  nonlinear  optics  where  the  index  of 
refraction  depends  on  intensity.  It  is  often  called  the  non¬ 
linear  Schrodinger ' s  equation.  The  PE  method  has  also  been 
applied,  in  more  recent  years,  to  the  field  of  fiber  optics 
and  plasma  physics. 

Beam  propagation  in  random  media  also  lends  itself  well 
to  the  PE  method.  Regardless  of  the  source  of  the  beam  -- 
radio,  acoustic,  optical  —  the  problem  is  analogous  to  the 
quantum  mechanics  problem  of  motion  in  a  random  potential. 
The  radar  application  of  this  problem  in  a  randomly  fluctu¬ 
ating  ionosphere  using  the  split-step  Fourier  algorithm, 
discussed  later  in  the  next  chapter,  is  covered  extensively 
by  Hardin  and  Tappert  (1974) . 
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The  most  extensive  use  of  the  PE  approximation  to  the 
elliptical  wave  equation  has  been  in  the  field  of  low-fre¬ 
quency  underwater  acoustic  propagation.  More  than  120  arti¬ 
cles  and  technical  reports  regarding  the  PE  method  and  its 
applications  to  underwater  acoustics  were  published  in  the 
last  15  years.  Numerous  marching  algorithms  and  computer 
applications  have  been  constructed  with  varying  degrees  of 
success.  The  earliest  of  these  were  based  on  the  split- step 
Fourier  algorithm  (Hardin  and  Tappert,  1974)  .  In  general, 
these  models  accept  input  defining  the  sound  speed,  volume 
loss  profiles,  and  depth  contours,  and  produce  output  of  the 
acoustic  field  and  associated  transmission  loss  (TL)  curves. 
These  results  are  generally  in  excellent  agreement  with 
experimental  measurements  for  deep  ocean  environments.  The 
theory  behind  this  method  and  subsequent  attempts  to  improve 
it  are  covered  in  detail  in  the  next  chapter.  Since  its  ini¬ 
tial  use,  the  PE  method  has  been  extended  to  include  higher 
acoustic  frequencies,  random  internal-wave  fluctuations  in 
the  index  of  refraction,  and  shallow-water  environments  with 
mixed  levels  of  success. 
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C.  ADVANTAGES  AND  LIMITATIONS  OF  THE  PE  METHOD 

Long-range,  low-frequency  sound  propagation  is  gener¬ 
ally  dominated  by  rays  having  small  grazing  angles  since 
rays  propagating  at  steep  angles  are  greatly  attenuated  due 
to  penetration  and  absorption  in  the  seabed.  Tappert  intro¬ 
duced  the  PE  method,  which  decomposes  an  elliptical  wave 
equation  into  two  equations  through  the  separation  of  outgo¬ 
ing  and  incoming  propagating  fields.  Tappert' s  implementa¬ 
tion  of  this  method  was  the  first  applied  to  underwater 
acoustics  and  is  referred  to  as  the  standard  parabolic  equa¬ 
tion  (SPE) .  The  SPE  implementation  and  many  of  its  more 
recent  incarnations  share  a  number  of  inherent  advantages 
and  limitations . 

The  advantages  of  the  PE  approximation  can  be  divided 
into  three  categories  (McDaniel,  1975)  .  First,  in  the  long- 
range,  low-frequency  environment  outlined  above,  the  govern¬ 
ing  equation  can  be  easily  expressed  by  a  parabolic  equa¬ 
tion,  which  is  easier  to  solve  than  either  the  elliptical  or 
the  hyperbolic  types.  The  PE  method,  therefore,  provides  a 
relatively  quick  solution  to  long-range,  low-frequency, 
range- dependent  problems.  Second,  in  solving  the  Helmholtz 
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equation,  the  problem  is  posed  as  a  pure  boundary  value  prob¬ 
lem.  In  practice,  the  determination  of  the  vertical  bound¬ 
ary  condition  is  often  difficult  and  imprecise.  The  PE, 
however,  poses  the  problem  as  one  of  an  initial-boundary 
problem,  therefore  avoiding  the  difficulty  of  the  vertical 
far-end  wall  boundary  condition.  Finally,  the  PE  solution 
is  generally  cheaper  to  obtain  in  both  memory  and  run-time 
than  solving  the  elliptical  equation. 

The  limitations  of  the  PE  method  can  be  placed  in  one  of 
two  broad  types.  The  first  of  these  results  from  the  mathe¬ 
matical  formulation  of  the  PE  itself  which  does  not  allow  the 
inclusion  of  certain  physical  phenomenon.  The  PE  is  only 
valid  in  the  far-field  and  therefore  cannot  be  used  for  a 
complete  analysis  of  near-field  effects.  The  index  of 
refraction  is  assumed  to  be  slowly  varying  in  range.  Large 
and  or  abrupt  changes  in  the  index  of  refraction  can  limit 
the  accuracy  and  validity  of  results.  The  PE  uses  the  one¬ 
way  wave  equation  and  wave  propagation  from  the  reverse 
direction  is  ignored.  Backscattering  effects  are,  there¬ 
fore,  neglected.  Treatment  of  the  so-called  "square  root 
operator"  defines  the  size  of  the  propagation  angle  which  in 
turn  effects  the  validity  of  the  PE. 
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The  second  type  of  limitation  is  defined  by  the  method 
of  solution  to  the  PE.  Many  of  the  limitations  of  this  type 
are  specific  to  the  solution  algorithm  being  implemented. 

The  fast  Fourier  transform  works  well  only  for  equations 
with  constant  coefficients,  and  therefore  must  be  used  with 
caution.  If  the  boundary  is  rigid,  the  algorithm  cannot 
treat  it  exactly.  Finite  difference  techniques  require  reg¬ 
ular  grid  partitions.  Using  the  horizontal-interface  PE  to 
handle  irregular  interfaces,  unless  the  sloping  angle  is 
small  or  density  change  is  small,  may  not  produce  satisfac¬ 
tory  results.  In  theory,  there  are  no  limitations  on  fre¬ 
quency.  In  practice,  however,  using  finite  difference 
techniques  for  high-frequency  problems  results  in  intolera¬ 
ble  and  expensive  computational  times . 

The  two  types  of  limitations  result  in  three  basic 
errors.  The  first  error  is  due  to  the  initial  PE  approxima¬ 
tion.  Findings  by  Smith  and  Smith  (1997)  and  others  show 
that  in  the  SPE  a  normal  mode  will  be  propagated  with  the 
correct  amplitude  and  mode  shape  but  with  an  error  in  phase 
velocity.  Some  higher  order  PE  approximations  improve  the 
accuracy  of  the  phase  velocity,  but  may  not  properly  propa¬ 
gate  the  mode  amplitudes.  The  error  in  phase  velocity  is 
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addressed  in  a  number  of  PE  approximation  methods  in  the  way- 
in  which  the  square  root  operator  is  expressed.  The  second 
error  is  introduced  by  the  selection  of  the  range  step  size. 
The  error  is  typically  found  in  the  second  or  third  order  of 
the  range  step  increment  (Jensen,  et .  al.,  1994).  If  the 
error  is  defined  as  En<  then  for  the  case  of  several  modes 
propagating,  the  error  produced  in  a  single  range  step  size 
would  be  E  =  ^En  .  The  third  type  of  error  is  due  to  trun- 

n 

eating  the  field.  The  initial  approximation  of  the  PE  poses 
basic  limitations  as  outlined  above.  The  errors  introduced 
by  truncating  the  field  are  rooted  mainly  in  computational 
efficiency  in  which  the  numerical  results  are  obtained.  A 
detailed  analysis  of  the  effects  of  range  step  choice  were 
investigated  by  Smith  (2000) . 

D.  THESIS  SUMMARY 

The  objective  of  this  thesis  is  to  explore  the  accuracy 
and  validity  of  various  approaches  to  the  PE  approximation 
method  as  it  applies  to  underwater  acoustics,  in  particular, 
how  each  of  these  methods  or  implementations  affects  the 
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errors  outlined  above.  The  resulting  accuracy  and  validity 
of  the  model  output  will  be  analyzed. 

Chapter  II,  Theoretical  Background,  reviews  the  theory 
relevant  to  the  material  in  the  subsequent  chapters .  Over¬ 
views  are  provided  of  the  parabolic  equation  approximation 
used  for  acoustic  propagation  modeling,  the  standard  para¬ 
bolic  equation  (SPE)  method,  the  wide  angle  parabolic  equa¬ 
tion  method  (WAPE) ,  and  the  reference  sound  speed  (cD) 
insensitive  parabolic  equation  method  (COIPE)  implementa¬ 
tions.  Additionally,  a  brief  overview  of  the  method  of  nor¬ 
mal  modes  is  discussed  with  emphasis  on  the  modal 
decomposition  of  acoustic  fields  generated  by  PE  models. 

Chapter  III,  Numerical  Implementation,  outlines  the 
numerical  modeling  techniques  used  in  implementing  the  cQ- 
insensitive  and  double-precision  WAPE  models. 

Chapter  IV,  Numerical  Results,  explores  the  results 
obtained  from  each  of  the  implementations  examined.  These 
include  SPE,  WAPE  (single  and  double  precision) ,  and  COIPE. 
Each  model  is  examined  for  validity  and  accuracy  against  a 
number  of  standard  test  cases.  The  models  were  also  analyzed 
for  their  ability  to  reproduce  the  standard  normal  mode 
basis  set. 
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Chapter  V,  Conclusions  and  Summary,  presents  a 
tion  of  conclusions  and  identifies  areas  of  further 
gation  and  research. 


suinma- 

investi- 
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II.  THEORETICAL  BACKGROUND 


This  chapter  outlines  the  theoretical  background  sup¬ 
porting  the  development  of  the  standard,  wide-angle,  and  c0- 
insensitive  implementations  of  the  parabolic  equation  (PE) 
approximation  to  the  acoustic  wave  equation.  The  bulk  of 
this  theoretical  development  is  directed  toward  its  applica¬ 
tion  in  the  Monterey-Miami  PE  (MMPE)  implementation  of  the 
wide-angle  parabolic  approximation.  The  MMPE  implementation 
is  emphasized  here  since  it  has  been  subjected  to  extensive 
analysis  and  numerous  numerical  improvements  over  its  life 
span.  It  also  serves  as  the  framework  of  the  COIPE  model 
discussed  and  analyzed  later.  A  summary  of  the  method  of 
normal  modes  with  an  emphasis  on  the  modal  decomposition  of 
acoustic  fields  generated  by  PE  models  is  also  included. 

The  entry  point  in  deriving  the  PE  is  the  Helmholtz 
equation,  in  cylindrical  coordinates, 


2  2 

I^PfrM  +  IiP_  +  ^P+k2n2D  _  A 
raiVaJ  r?a<p2  az2  0  p "  ' 


(2.1) 


where  p(r,z)  is  the  acoustic  pressure,  k0  =  od/c0  is  the  ref¬ 
erence  wavenumber,  n(r,  z,  <p)  =  c0/c(r,  z,  <p)  is  the  acoustic  index 


of  refraction,  c0  is  the  reference  sound  speed,  and  c(r,  z,  cp) 
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is  the  acoustic  sound  speed.  All  of  the  environmental  char¬ 
acteristics  are  represented  in  n(r,  z,  <p) .  The  treatment  of 
density  contrasts  will  not  be  presented  here,  but  is 
included  in  the  MMPE  model  by  an  effective  index  of  refrac¬ 
tion  term  (Smith,  2000) .  A  time-harmonic  component  of  the 
acoustic  field  then  has  the  form 

P(r,  z,  (p,  rot)  =  p(r,  z,  (p)e-ia)t .  (2.2) 

A.  OPERATOR  NOTATION 


In  order  to  develop  the  various  parabolic  approxima¬ 
tions  to  the  Helmholtz  equation,  it  is  useful  to  introduce  an 
operator  notation.  Let 


P 


d_ 

°P  dr 


(2.3) 


and 


Qod  =  (p  +  e  +  u  +  1)172. 


^op 


(2.4) 


.  2  i  1  d  2  .  1  dd 

where  £  =  nr  -  1  ,  p  =  — >  and  u  = 


k$dz2 


kh2d(p2  ’ 


To  account  for 


the  dominate  cylindrical  spreading  and  to  simplify  the  form 
of  the  Helmholtz  equation,  the  acoustic  pressure  may  be 

defined  as  p(r,  z)  =  -^q(r,  z),  which  yields 

Jr 


14 


(P»p  +  ko2Qop)q  =  0  •  (2-5) 

Proper  factorization  of  this  equation  is  accomplished  by¬ 


defining  q  =  Q~*/2u  ( Tapper t,  1977)  ,  such  that 

(Pop  +  ikoQop)(Pop  -  ikoQop)u  +  iko[pop’  Qop]u  =  0  •  (2.6) 

The  commutator  [Pop,Qop]  is  exactly  zero  for  layered  media, 
and  is  assumed  negligible.  The  remaining  factors  of  Eq. 

(2.6)  define  the  wave  equations  for  the  incoming  and  outgo¬ 
ing  fields,  the  latter  of  which  then  satisfies 

Popu  =  ikoQopu  (2-7) 

or 

-ik0^  =  Q°PU  •  (2-8) 

In  cases  where  backscattered  energy  may  be  neglected,  Eq. 
(2.8)  provides  the  full  description  of  the  forward  propagat¬ 
ing  acoustic  energy. 

B.  THE  SPLIT- STEP  FOURIER  ALGORITHM 

The  most  common  methods  of  computing  PE  solutions  are 
(1)  the  finite  element  method  (Collins,  1994),  (2)  the 

implicit  finite  difference  method  (Lee,  et.  al.,  1981),  and 
(3)  the  split-step  Fourier  (PE/SSF)  method  (Hardin  and  Tap- 
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pert,  1973) .  Of  these  three  methods,  the  PE/SSF  method  has 
the  distinct  advantage  of  speed.  This  is  accomplished  by- 
decomposing  the  acoustic  field  into  a  slowly  modulating 
envelope  function  and  a  phase  term  which  oscillates  with  the 

acoustic  frequency.  Defining  u  =  'Felk°r  and  substituting  into 
Eq.  (2.8)  yields 

=  -  *0*  +  *oQop'*'  =  -ikoHop't' .  (2.9) 


where  Hop  =  l-Qop  is  a  Hamiltonian- like  operator  which 

defines  the  evolution  of  the  PE  field  function  in  range. 

In  order  to  compute  the  solution  to  the  acoustic  field, 
a  marching  algorithm  must  be  developed  of  the  form 

4/(r  +  Ar)  =  <I>(r)'F(r)  .  (2.10) 

The  propagator,  <I>(r) ,  is  designed  to  propagate  the  solution 
out  in  range.  The  PE/SSF  method  accomplishes  this  via 

<I>  ( r ) «  e-ikoH°p(r)Ar  (2.11) 

where 


i  r<r*4r) 

Hop(r)  =  Hop(r’)dr\  (2.12) 

r 

In  practice,  the  Hamiltonian  is  usually  considered  constant 
over  a  small  range-step  such  that  H0p(r)  =  Hop(r)  . 
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The  operator  Hop  (and  Qop)  is  not  a  scalar  operator  but 


a  combination  of  scalar  and  differential  operators.  Each 
individual  operator  within  Hop  can,  however,  be  applied  by  a 
simple  scalar  multiplication  in  the  proper  domain.  There¬ 
fore,  each  of  the  terms  of  Hop  must  be  separated  in  order  to 
comply  with  the  PE/SSF  algorithm.  This  necessitates  the 


approximation  of  the  square-root  operator,  specifically 

(2.13) 


2  2  2  2 
H°p(^2’  n(r’  Z))  “  T°p(a^2)  +  U°P(n(r’  z))  +  vop(ria^) 


It  is  this  approximation  that  forms  the  basis  of  numerous 
higher-order  PE  model  implementations. 

If  the  uncoupled  azimuth  approximation  is  employed 
then  Vop  =  0  identically.  The  Uop  is  a  multiplication  opera¬ 
tor  in  z-space  and  is  therefore  a  diagonal  matrix.  The  oper¬ 
ator  Top  is  not  diagonal  in  z-space,  resulting  in  different 
eigenfunctions  being  coupled.  The  Top  is,  however,  diagonal 
in  vertical  wavenumber  space.  It  is  advantageous  then  to 
treat  each  operator  separately,  one  in  kz-space,  and  the 
other  in  z-space.  Employing  the  Baker-Campbell-Hausdorf f 
expansion  (Bellman,  1964), 

gA  +  B  _  gAgBg[A,  B]  +  [A,  [A,  B]]  +  [B,  [B,  A]]  +  ...  ,  (2.14) 
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Top  and  Uop  are  both 


where  A  =  -ik0ArTop  and  B  =  -ik0ArUop  . 
small,  therefore  it  is  assumed  their  products  of  second 
order  are  negligible.  Using  Eq.  (2.14),  <I>(r)  becomes 


<D(r)  =  e 


_  e  lk°2  U°P(r  +  Ar)  -ik0ArT  lko  2  UoP(r) 


(2.15) 


Third  order  accuracy  results  from  this  centered  step  scheme 
(Jensen,  et.  al . ,  1994).  Due  to  the  formulation  of  the  prop¬ 
agator  there  are  no  intrinsic  losses  due  to  the  numerical 
scheme  and  energy  conservation  is  retained  (Tappert,  1995) . 

In  summary,  the  PE/SSF  algorithm  accepts  the  PE  field 
function  'F  specified  at  some  range  r  in  the  z -domain.  This 

-ik0^-Uop(r) 

is  followed  by  a  multiplication  of  e  2  ,  the  z-space 

operator  defined  at  the  beginning  of  the  range -step.  This 


product  is  then  transformed  to  the  kz-domain  and  multiplied 


by  e  lk°ArTop(kz) ,  the  k2-space  operator.  Another  transformation 
to  the  z-domain  and  multiplication  by  the  z-space  operator 

-ik0^Uop(r  + Ar) 

e  2  ,  defined  at  the  end  of  the  range-step,  produces 

the  final  resultant  field  function  at  r+Ar  in  the  z-domain. 
Note  that  the  FFT  convention  used  here  and  in  the  numerical 
implementations  is 
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(2.16) 


Y(z)  =  FFT(xF(kz))  , 

T(k7)  =  IFFT('F(z))  * 

Wrapping  these  steps  into  a  single  equation  results  in 


'P(r  +  Ar,z)  = 


(2.17) 


-ik0yUop(r  +  Ar,z) 


x  FFtJ e-ikoArT»p(kz)  x  IFFlfe  lk° 2  U°p(r’ Z)  x 


W(T,  2)] 


C.  STANDARD  PE  DEVELOPMENT 


The  lowest  order  approximation  to  the  Hamiltonian  oper¬ 
ator  is  commonly  referred  to  as  the  standard  parabolic  equa¬ 
tion  (SPE) .  This  approximation  is  obtained  by  assuming  the 
operators  are  small. 


£ « 1  and  ja  «  1  .  (2.18) 

The  first  of  these  is  merely  a  recognition  that  sound  speeds 
vary  by  less  than  2%  in  most  ocean  environments,  so  n(r,  z)~\ 
everywhere.  The  second  assumption  can  be  seen  to  imply  that 


vertical  angles 


^k2 

,k0 


=  sin20 


are  small.  This  limits  the  SPE 


accuracy  to  small  angles. 

With  these  approximations,  the  SPE  is  obtained  by  a 
binomial  expansion  of  the  square  root  operator, 
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(2.19) 


H 


spe 


The  separated  operators  then  have  the  form 

Tspe(d^  =  ~2kfizi  (2,20) 

and 

Uspe(r,  z)  =  -^(n2(r,  z)  -  1 )  .  (2.21) 

Note  that  in  the  context  of  Hspe  acting  as  a  Hamiltonian 
operator,  the  operators  Tspe  and  Uspe  correspond  to  kinetic 
and  potential  energy  operators,  respectively.  In  the  verti¬ 
cal  wavenumber  domain, 

lkz 

Tspe(kz)  =  .  (2.22) 

2k0 


D.  WIDE-ANGLE  PE  DEVELOPMENT 


A  higher  order  approximation  to  the  operators  was 
introduced  by  Thompson  and  Chapman  (1983)  and  is  often 
referred  to  as  the  widerangle  PE  (WAPE)  approximation.  The 
operators  of  the  WAPE  are  defined  by 


LwaP‘W, 


=  1  - 


21 


1  + 


±JL 

kfidz2 


1/2 


(2.23) 
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and 


Uwape(r’z)  =  ~(n(r,  z)  -  1)  .  (2.24) 

In  wavenumber  space,  the  WAPE  kinetic  energy  operator  takes 
the  form  of 


Twape(kz) 


1- 


2-,  1/2 


(2.25) 


Note  that  energy  propagating  at  angles  beyond  vertical, 
kz>k0,  is  evanescent  since 


Twape(kz  >  ko) 


(2.26) 


The  wide-angle  PE,  and  other  higher-order  approxima¬ 
tions  to  the  SPE,  was  originally  developed  in  an  effort  to 
extend  the  region  of  validity  of  the  SPE.  While  the  wide- 
angle  PE  succeeded  in  this  goal,  it  introduced  a  number  of 
other  issues.  These  include  a  sensitivity  to  the  selection 
of  the  reference  sound  speed  and  an  inability  to  reproduce 
the  standard  normal  mode  basis  set.  This  said,  however,  the 
wide-angle  approximation  has  been  very  successful  in  produc¬ 
ing  useful  and  valid  results  in  many  situations.  The  follow¬ 
ing  section  addresses  an  attempt  to  mitigate  the  reference 
sound  speed  sensitivity  experienced  by  the  WAPE  approxima¬ 
tion. 
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E.  COIPE  DEVELOPMENT 


WAPE  and  other  higher-order  models  are  not  exact  repre¬ 
sentations  and  may  under  certain  circumstances  produce 
results  worse  than  the  SPE  which  they  are  meant  to  improve 
upon  (Porter  and  Jensen,  1993)  .  In  these  cases,  it  is  often 
found  that  the  error  results  from  errors  in  the  underlying 
normal  mode  basis  set,  their  associated  phase  speeds,  and  a 
sensitivity  to  the  choice  of  reference  sound  speed,  c0.  The 
choice  of  c0  is  the  one  ambiguous  feature  of  all  PE  models. 
Due  to  this  ambiguity,  it  is  desirable  to  develop  a  model 
that  is  insensitive  to  the  input  reference  sound  speed. 

Then,  any  significant  deviations  in  the  choice  of  c0  would 
not  adversely  affect  the  computed  results.  In  an  attempt  to 
develop  a  wide  angle  implementation  that  is  insensitive  to 
c0,  Tappert  (1991)  developed  the  c0-insensitive  PE  (COIPE) 


approximation . 

The  COIPE  approximation  based  on  the  WAPE  approximation 
with  the  operator  replacement 


_1_  02  1  0  _j.  .10 

kfidz1  ko(r)0zn  f’  Z'ko(r)0z  ' 


(2.27) 
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where  k0(r)  is  now  range-dependent  due  to  the  range- dependent 


nature  of  the  reference  sound  speed,  c0(r)  .  Introducing  the 
operator 


p(r)  =  - 


i  a 

k0(r)3z  ' 


(2.28) 


the  c0-insensitive  approximation  to  the  Hamiltonian  operator 
becomes 


HIns(r,z)  =  ( 1  -  n(r,  z))  +  ( 1  -  a/i  -  p(r)n_1  (r,  z)p(r))  .  (2.29) 

In  this  form  the  spherical  wavefront  property  of  the  TC 
approximation  is  maintained  for  steep  angles.  Furthermore, 
this  approximation  is  not  sensitive  to  the  choice  of  c0(  and 
reduces  to  a  c0-independent  approximation  for  small  angles 
(Tappert,  1977) .  The  resulting  wave  equation  approximation 
is 


<)¥ 

5jT  =  -^0Hlns(r)T  . 


(2.30) 


To  implement  this  approach,  a  transform  to  a  new  depth 
coordinate  is  necessary.  The  transformation  relationships 


are 


z 

z  =  J Vn(z,  r’)dr’ , 

o 


(2.31) 
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n(z,  r)  =  n(z,  r) , 


(2.32) 


and 


T(z,r)  =  ir1/4(z,r)T(z,r).  (2.33) 
The  corresponding  Hamiltonian  operator  then  becomes 

H(r)  =  (l-n(z,r))  +  (l-Vl-p2(r))  .  (2.34) 
Note  that  in  the  tilde-domain  this  has  the  same  form  as  the 
TC  operator.  The  resulting  form  of  the  parabolic  wave  equa¬ 
tion  in  the  tilde  domain  is 


=  ik0Hins(r)'P  .  (2.35) 

The  corresponding  kinetic  and  potential  energy  operators  of 
the  COIPE  in  the  tilde-domain  may  then  be  represented  by 


Tins(r,  k~) 


li__ Ljl 

q  tf(r)dz2 


(2.36) 


and 


Uins(r,  z)  =  1  -  n(r,  z)  .  (2.37) 

The  COIPE  requires,  in  general,  the  computation  of  c0 
for  each  range  step.  To  compute  c0,  the  condition  that  the 

water  depth  be  unchanged  in  the  z  coordinate  must  be  satis¬ 
fied.  This  is  accomplished  via 
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where  h(r)  is  the  water  depth  at  range  r.  Under  this  condi¬ 
tion,  the  reference  sound  speed  is 

/  h(r) 

C°(r)  _  h(r)  I  7^7) 

^  o 

The  University  of  Miami  PE/SSF  acoustic  model  (UMPE) 
has  been  tested  using  the  technique  outlined  above  with 
excellent  results  (Tappert,  1995) .  A  numerical  implementa¬ 
tion  of  the  COIPE  within  the  existing  framework  of  the  MMPE 
model  is  covered  in  the  following  chapter.  Comparative 
analysis  of  this  approach  with  the  SPE,  MMPE,  and  selected 
normal  mode  test  cases  is  covered  in  Chapter  IV. 

F.  NORMAL  MODE  FUNCTIONS  AND  THEIR  DECOMPOSITION  FOR  THE 
PE  APPROXIMATION 


One  of  the  difficulties  encountered  with  the  WAPE 
approximation  is  that  it  does  not  decompose  into  the  stan¬ 
dard  normal  mode  basis  set  as  defined  by  normal  mode  theory 
and  the  SPE.  Thomson  (1993)  looked  at  this  problem  for  the 
SPE  approximation.  He  showed  that  the  field  p  that  satisfies 


the  acoustic  wave  equation  exactly  relates  to  the  field 
satisfying  the  SPE  approximation.  The  wavenumber  of  W  and 
the  modal  amplitudes  may  be  determined  exactly.  The  corre¬ 
sponding  normal  mode  amplitudes  and  wavenumbers  are  then 
obtained.  The  WAPE,  however,  has  modal  eigenfunctions  and 
eigenvalues  of  distinct  from  the  acoustic  wave  equation. 

Smith  and  Smith  (1997)  showed  that  the  modal  decomposi¬ 
tion  of  solutions  based  on  the  WAPE  generated  erroneously 
range -dependent  mode  amplitudes  in  a  range- independent  envi¬ 
ronment.  This  was  due  to  the  mismatch  between  the  standard 
normal  mode  basis  set  (based  on  the  Helmholtz  equation)  and 
the  proper  basis  set  of  the  WAPE.  Furthermore,  they  were 
able  to  develop  an  approximate  numerical  scheme  for  comput¬ 
ing  the  proper  WAPE  basis  set.  This  resulted  in  the  proper 
range -independent  character  of  the  mode  amplitudes.  Addi¬ 
tionally,  they  suggested  that  the  higher-order  COIPE  version 
of  the  WAPE  might  overcome  such  issues  and  decompose  prop¬ 
erly  into  the  standard  normal  mode  basis  set.  This  will  be 
one  of  the  primary  tests  of  the  improved  accuracy  of  the 
COIPE  over  the  WAPE. 

Normal  modes  are  obtained  by  using  the  method  of  sepa¬ 


ration  of  variables.  As  presented  here,  it  assumes  the  ocean 
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is  horizontally  stratified,  is  of  constant  water  density,  is 
azimuthally  symmetric,  and  is  range- independent .  As  was 
presented  for  the  PE,  the  time  harmonic  acoustic  field 

p(r, z)e_10)t  at  r>0  due  to  a  point  source  at  z  =  zs  and  r  =  0  sat¬ 
isfies  the  2D  Helmholtz  equation  (no  azimuthal  variations 
possible) 

rKr5?)+|#  +  k»2n2(z)P  =  0'  ,2-40) 


For  solutions  of  the  form  p(r,  z)  =  —  <|>(r)v(z) ,  the  depth- 

Jr 


dependent  modal  equation 

.2 

-^-vm(z)  +  (k2n2(z)-K2)vm(z)  =0  (2.41) 

dzz 

is  obtained.  Here,  ,  the  square  of  the  horizontal  wave- 
number  for  mode  m,  is  the  separation  constant  (eigenvalue) 
and  vm(z)  is  the  specific  mode  function  (eigenfunction) 

associated  with  the  separation  constant.  Given  an  arbitrary 
pressure  field,  the  sum  of  the  normal  modes  is  required  and 
is  expressed  as 


p(r, z)  =  -j  £  <t>m(r)vm(z) ,  (2.42) 

Vr 

m  =  1 
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where  the  cylindrical  spreading  factor  has  been  explicitly 
included.  The  factor  <t>m(r)  can  easily  be  shown  to  have  the 

form  of  a  Hankel  function  scaled  by  mode  excitation  values 
depending  on  the  source  strength  and  depth  in  the  waveguide. 
In  the  far- field,  the  only  range- dependence  of  this  function 
is  in  the  phase  (having  already  accounted  for  cylindrical 
spreading).  Thus,  a  generalized  mode  amplitude,  Am(r) ,  may 
be  defined  according  to 

p(r,  z)  =  Am(r)vm(z)  .  (2.43) 

m 

The  normalized  eigenfunctions  form  a  complete,  orthogo¬ 
nal  basis  set,  defined  by 

D 

JVn(z)Vm(z)dz  =  5mn  •  (2  •  44 ) 

0 

Thus,  by  using  the  orthogonality  relation  and  assuming 
normalized  modes,  the  generalized  modal  amplitudes  can  be 
extracted  from  the  computed  pressure  field  by  forming  an 
inner-product 

D 

A  /  X  r  fp(r,  z)vm(z) 

Am(r)  =  VrJ - - dz.  (2.45) 

0 
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This  technique  is  used  in  the  analysis  and  comparisons  of  the 
normal  mode  amplitudes  of  the  SPE,  WAPE,  and  COIPE  models  in 
Chapter  IV. 

Finally,  the  corresponding  depth  separated  equations  of 
the  various  parabolic  approximations  introduced  should  be 
examined.  First,  note  that  the  normal  mode  equation,  Eq. 
(2.41),  may  be  written  in  terms  of  the  operators  |i  and  e 
(defined  in  Eq.  (2.4))  as 

(l*  +  e)vm  =  ("p7  “  l)vm  •  <2.46) 

This  indicates  that  the  basis  set,  vm  ,  are  eigenfunctions  of 


the  combined  operator  (n  +  £)  with  eigenfunctions  (K2m/k20-l). 

The  PE  approximation  has  a  depth  separated  equation  of 
the  form 


K 


H  y  =  i— y 

op/wn  V  ' 
K0 


(2.47) 


where  the  functions  xm  are  the  eigenfunctions  of  the  corre¬ 
sponding  PE  approximation.  In  order  for  the  eigenfunctions 
Xm  and  vm  to  be  the  same  (up  to  some  normalization  factors)  , 

the  operator  Hop  must  be  some  rational  function  of  the  oper¬ 


ator  (n  +  e).  For  example,  a  simple  polynomial  series  with 
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terms  (n  +  £)n,  where  n  is  some  integer,  will  have  the  same 


basis  set  vm . 

The  standard  parabolic  equation  is  a  simple  first  order 
binomial  expansion  of  the  Hamiltonian  operator,  i.e. 

HsPe  =  +  e)  •  (2.48) 

This  obviously  satisfies  the  above  criterion,  therefore  the 
SPE  approximation  has  the  same  normal  mode  basis  set  as  the 
fundamental  Helmholtz  equation. 

The  wide-angle  PE,  on  the  other  hand,  may  be  written  in 
the  form 

Hwape  =  2-Vl+H-Vl+e.  (2.49) 

Performing  a  power  series  expansion  on  this  operator  to  sec¬ 
ond  order  yields 

Hwape“-|(^+e)  +  §(^2  +  e2)  •  (2.50) 

This  shows  that  the  WAPE  approximation  is  only  first  order 
accurate  in  the  operator  (|X  +  s)  but  has  errors  in  the  second 
order  terms.  It  is  this  weakness  of  the  WAPE  which  generates 
some  of  the  associated  problems  with  this  approximation. 

The  c0-insensitive  PE  approximation,  however,  may  be 
written  as 
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(2.51) 


A  power  series  expansion  of  this  operator  to  third  order 
yields 

Hi„s“-|(l»  +  e)  +  |((i  +  e)2-i(|i  +  e)3  +  l^.  (2.52) 

Thus,  the  COIPE  is  accurate  to  second  order  in  (ji  +  e)  with 
leading  error  in  the  third  order  terms.  This  was  first  shown 
by  Tappert  and  Brown  (1996)  .  It  is  this  feature  of  the  COIPE 
which  makes  it  an  attractive  alternative  to  the  WAPE  approx¬ 
imation.  As  will  be  shown,  this  also  generates  solutions 
which  may  be  decomposed  into  the  standard  normal  mode  basis 
set  without  significant  error. 
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III.  NUMERICAL  IMPLEMENTATION 


This  chapter  outlines  the  numerical  implementation  of 
the  PE  models  used  in  the  analysis  of  the  c0-insensitive  PE 
approximation  (COIPE)  as  well  as  a  double-precision  version 
of  the  MMPE.  The  MMPE  (Smith,  2000)  is  an  existing  well-doc¬ 
umented  version  of  the  WAPE  approximation.  In  the  cases 
where  a  SPE  solution  was  required  for  comparison  purposes, 

the  operators  Uwape  and  Twape  were  replaced  with  Uspe  and  Tspe 

within  the  framework  of  the  MMPE  model.  As  was  discussed  in 
the  previous  chapter,  the  COIPE  implementation  strives  to 
eliminate  the  sensitivity  to  the  chosen  reference  sound 
speed  suffered  by  the  WAPE  approximation.  As  will  be  shown 
in  the  next  chapter,  a  small  change  in  the  specified  refer¬ 
ence  sound  speed  can  result  in  a  dramatic  change  in  the  pre¬ 
dicted  acoustic  field.  This  error  often  increases  with 
range  and  for  complex,  highly  range- dependent  environments 
can  quickly  result  in  the  produced  data  being  unusable.  The 
numerical  implementation  described  here  was  chosen  such  that 
it  could  be  easily  integrated  into  the  existing  MMPE  model  to 
provide  the  most  reliable  comparison  of  results.  The 
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results  of  these  comparisons  are  discussed  in  the  following 
chapter . 

A.  ENVIRONMENTAL  GRID  SIZES 

The  MMPE,  as  well  as  the  COIPE,  implementation  dis¬ 
cussed  here  uses  a  discretized  grid  to  describe  the  environ¬ 
ment.  The  grid  step-sizes  are  user  definable  and  can 
dramatically  affect  the  solution  obtained  by  the  model.  The 
grid  sizing  method  described  here  is  that  used  by  the  MMPE 
model  as  well  as  the  COIPE  implementation  of  the  MMPE.  The 
environment  is  discretized  by  a  mesh  of  size  (Ar,  Az).  This 
results  in  the  field  and  propagator  functions  being  dis¬ 
cretized  in  depth  with  array  length  N,  where  N  is  the  corre¬ 
sponding  length  of  the  FFT  used  in  the  PE/SSF  algorithm.  In 
other  words, 

'F(r,z)=*'F(ri,zn),  (3.1) 

where 

r*  =  (i-l)A,  i  =  1  ...imax  ,  (3.2) 

and 
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zn  = 


n.  ,  n 

n--\Az,  n  =  1»2 


-(n-„4)Az-  n  =  f  +  1> 


(3.3) 


N 


The  depth  mesh  is  defined  such  that  grid  points  lie  on  frac¬ 
tional  values  of  Az .  This  convention  results  in  the  avoid¬ 
ance  of  carrying  the  zero-pressure  value  at  z  =  0  throughout 
the  calculations.  It  should  be  noted  that  because  of  using 
the  full  Fourier  transform,  half  of  the  depth  mesh  values 
define  an  image  ocean  for  negative  depths.  This  has  the 
added  benefit  of  enforcing  the  surface  boundary  condition, 
defined  in  the  next  section,  via  FFT  symmetry.  Analysis  by 
Tapper t  (1977)  and  Smith  (2000)  showed  that  by  considering  a 
lower  limit  on  the  allowable  angles  of  propagation,  a 
default  value  of  Az  may  be  defined  and  therefore  a  default 
transform  size  N.  This  is  possible  since  the  depth  mesh 


influences  the  wavenumber  increments  Akz  via  the  FFT.  For  a 


given  environment  with  a  relatively  small  angle  of  propaga¬ 
tion,  the  mesh  size  may  be  increased.  An  increase  in  the 
mesh  size  results  in  a  corresponding  decrease  in  the  compu¬ 
tational  time.  As  propagation  angles  increase,  the  mesh 
size  is  reduced.  The  most  accurate  solution,  therefore, 
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should  be  obtained  when  both  Az  and  Ar  are  on  the  order  of  a 
wavelength. 


B.  TILDE  TRANSFORMATION  AND  C0  GENERATION 


The  tilde  transformation  begins  by  defining  a  rescaled 
depth  variable  z  in  terms  of  the  true  depth  z 


z 


0 


If  it  is  required  that  the  bottom  depth  remains  fixed,  in 
other  words  that  Zb  =  zb  =  h  ,  then  the  range -dependent  bottom 
depth  is 

h(r)  h(r)  _  h(r) 

h(r)=  \  =  f  J^F)  f  -p!=dz’,  (3.5) 

J  J  A/c(z,  r)  N  u  J  7c(z’,  r) 

0  0  0 

where  c0  has  been  specifically  declared  as  range -dependent, 

due  either  to  h(r)  or  c(z,r).  Solving  for  c0(r). 


c0(r)  = 


1 


h(r) 

—  f - 

h<r)  j  Jc(z\  r) 
0 


dz’ 


-2 


(3.6) 


This  expression  defines  the  rang-dependent  reference  sound 
speed  to  be  used  in  the  parabolic  approximation.  Its  decla- 
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ration  then  eliminates  the  need  for  a  user-defined  value. 
Furthermore,  since  the  definition  is  based  on  an  evaluation 
of  the  local  environment  at  each  range,  it  will  preserve  rec¬ 
iprocity. 

If  c(z)  varies  linearly  with  depth, 

c(z)  =  c(zi_1)  +  b(z-zi_1)  =  a  +  bz,  (3.7) 

where  b  is  the  sound  speed  gradient, 


b  = 


c(z)-c(zi_1) 


(3.8) 


and 


a  =  c^i.^-bZj.j  , 

then  the  integral  may  be  written 


(3.9) 


f  -p=dz’  =  f  1  dz» . 
j  Vc(z')  J  ./n  +  hy’ 


Va  +  bz’ 


Z>-1 


Using  a  table  of  integrals,  this  becomes 


,  Ja  +  bz’ 


2  I _ 

=  jit  Va  +  bzj  —  J&  +  bZj_  j] 


-i-i 


(3.10) 


(3.11) 


Substituting  into  Eq.  (3.11)  for  a  and  b, 


J 


VcCzO 


:dz’  = 


2(zi-zi_1) 

(ci-ci-i) 


Vci  -  1 )  - 


2(zi~zi-l) 


(3.12) 
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Now  consider  the  sampling  of  the  environment  within  the 


model.  The  typical  definitions  are 


zk 


# 


k  = 


1,2,3... 


N 


$ 


(3.13) 


where  N  is  the  FFT  size.  Based  on  this  grid  sampling  of  the 
environment,  it  would  then  be  unusual  if  znbi  =  h(r)  .  Rather, 
the  bottom  lies  just  below  the  index 


(3.14) 


such  that  znbi<h(r).  To  perform  the  integrals  needed,  how¬ 
ever,  the  conditions 

zj  =  0  ,  (3.15) 

zNs  =  h(r)  / 

must  be  satisfied.  This  may  be  accomplished  by  first  trans¬ 
forming  to  a  Zj  space  (which  is  not  equally  spaced)  ,  which 

satisfies  Eq.  (3.15),  and  then  interpolating  Zj  to  z^  where 

Zk  =  (k-^Az.  (3.16) 

Part  of  the  difficulty  is  to  define  the  sound  speed  on 
the  non-grid  points  z  =  0  and  z  =  h(r)  .  Fortunately,  this  is 
essentially  already  accomplished  in  the  MMPE  model.  The 
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surface  sound  speed  is  required  as  an  input  to  the  MMPE,  and 


is  stored  in  the  sound  speed  array  by  the  environmental  prop¬ 
agator  subroutine.  The  sound  speed  at  the  bottom  depth, 
h(r),  is  also  already  estimated  at  the  grid  point  nbi .  While 
this  is  an  estimate,  it  provides  a  good  approximation  for 
cases  where  energy  doesn't  significantly  interact  with  the 
bottom.  This  results  in 

zj  =  0=>z  =  0=>c(0)  =  SS  ,  (3.17) 

zNs  =  h(r)=>z  =  h(r)  =>  c(h(r))  =  SSBW  , 

where  SSO  is  the  surface  sound  speed  input  into  MMPE  and  SSBW 

is  the  estimated  water  sound  speed  at  the  water/bottom 

interface . 

The  numerical  scheme  for  computing  c0  and  the  tilde 
transformation  can  now  be  defined.  Using  Eqs.  (3.6)  and 
(3.12) , 


c0(r)  =  (3.18) 


h2(r) 

zi 

i 

4  1 

N 

7? 

! 

N 

1 

H- i 

1 

h(1")  —  znbi 

l"2 

4 

7^  +  Vsso. 

*T 

_k  =  2'^k  +  VCk-l_ 

T 

[vssbw + VcnblJ 

*• 
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where  it  is  assumed  that  the  sound  speed  has  already  been 
interpolated  onto  the  gridded  mesh  zk  at  values  ck.  Simi¬ 
larly,  Eq.  (3.4)  with  Z\  =  0  becomes 


z2=2  Jc0,  (3.19) 

Zj  =  Zj_i  +2jro-=^ — j=  3...nbi  +  1  , 

V<Wck-i 

:  :  rr  ^ -  znbi 

zNs  =  zNs-  1  +  2  A/Cn 

VSSBW  + 

Note  that  Ns  =  nbi  +  2  since  the  points  z  =  0  and  z  =  h  have 
been  added.  As  a  result  of  this  transformation  and  the  above 
definitions , 

Cj  =  SS,  (3.20) 

CNs  =  SSBW, 
cj*l,  Ns  =  cj_i 

Given  these  relationships,  it  is  now  necessary  to 
interpolate  in  depth  the  tilde  domain  sound  speed  onto  the 
standard  grid.  Defining 

zk  =  (k-^Az,  (3.21) 

k=  1  2  ^  , 

IV  —  a,  2 

it  is  desired  that 

c(zk)  =  c(zk)  .  (3.22) 
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A  simple  linear  interpolation  from  c(zj)  to  c(zk)  could  be  per¬ 
formed.  This,  however,  would  not  preserve  the  linear  gradi¬ 
ent  nature  of  c(zk) ,  and  the  relation  defined  by  Eq.  (3.22) 
would  not  be  satisfied.  Instead,  assume 

zj_i<zk<Zj,  (3.23) 

and  then  define  the  fractional  depth  as 

fk  =  .  (3.24) 

zj  -  zj  - 1 

Tappert  (1995)  shows  that  the  inversion  of  the  mapping 
needed  produces 

tk  =  Cj _ i  +  gk (c,j  —  Cj - 1 ) ,  (3.25) 

where 

gk  =  r  +  {ft- (3.26) 

V Cj  +  VCj  -  1 


C.  ACOUSTIC  FIELD  GENERATION 


In  the  MMPE  model  the  source  amplitude  is  defined  rela¬ 
tive  to  a  small  finite  distance  from  the  source.  This  is 
done  to  avoid  the  singularity  at  range  r  =  0  in 


P 


fR0 

=  P°aT7 

N  r 


-'Feik°r . 


(3.27) 
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The 


Therefore,  the  source  is  defined  such  that  p(r  =  R0)  =  P0  . 

reference  range  is  defined  to  be  one  meter  to  remain  consis¬ 
tent  with  most  sonar  equations.  Using  this  convention,  the 
source  level  is  related  to  P0  by 

SL  =  201og(^  dB  re  Pr  R0  (3.28) 

where  Pr  =  1  pPa  at  the  reference  range  of  R0 - 


The  form  of  the  source  field  ^(r  =  0,  z)  must  now  be 
determined.  Rewriting  Eq.  (3.27)  as 

^(r,  z)elk°r  =  &>(r,  z)  ,  (3.29) 

F(WR0 

then 


xF(r  =  0,z)=  lim  =|-  E p(r,z), 

r— >0  PqVR0 


(3.30) 


where  in  the  vicinity  of  the  source  the  pressure  filed  takes 


the  form  of  the  spherical  Green's  function.  If  Ipl  =  — ^  is 

r  4ji 


required  at  R  =  R0 ,  then  let  a  =  P0R0  an<3  therefore 


p  =— eik°R 

p  4nR 


(3.31) 


where  R  =  Jr2  +  z2  .  The  source  is  represented  as  a  point 
source  at  (0,zs)  by 
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T(r  =  0, z)  =  a8(z - zs)  .  (3.32) 

Integrating  Eq.  (3.30)  over  all  depths  of  the  real  ocean 


results  in 


oo 


1 

4tiR 


e 


ik°rdz  . 


(3.33) 


0 

In  the  far-field,  r»z,  R  may  be  approximated  such  that 

R  =  Jr2  +  z2  =  r  +  z  =  z-zs  (3.34) 

thus  reducing  Eq.  (3.33)  to 


oo 


o 


and  leaving 


a  = 


(3.36) 


Performing  a  Fourier  transform  of  Eq.  (3.32)  with  the 
influence  of  the  image  source  included,  and  specifying  the 
source  in  the  kz-domain  yields 


'R(r  =  0, k)  =  aj  [8(z-zs)-8(z  +  zs)]e-ik*zdz  =  -2iasin(kzzs)  .  (3.37) 
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Equation  (3.37)  states  that  the  starting  field,  in  the  wave- 
number  representation,  has  constant  amplitude  modulated  by  a 
phase  due  to  the  interaction  of  the  source  and  its  image. 
This  representation  is  consistent  with  an  omnidirectional 
source  that  supplies  equal  energy  to  all  wavenumbers . 

Rather  than  setting  the  amplitude  of  this  function  to 
unity  to  ensure  equal  population  of  all  directions,  a  smooth 
taper  is  used  at  high  absolute  wavenumber  values  to  limit  the 
angular  dimension  of  the  source  and  to  reduce  sidelobe 
influences.  This  is  necessary  due  to  the  wide-angle  approx¬ 
imation  being  valid  to  only  approximately  40°  ,  and  due  to  the 
restriction  on  how  large  the  finite  FFT  will  allow  kz/kQ . 

The  wide-angle  source  is  modified  to  produce  an  accurate 
solution  in  the  far- field  (Thomson  and  Bohun,  1988)  by 

(  k2V,/4 

F(kz)  =  l--§  ,(|kj<k0).  (3.38) 

v  ko  J 

The  limits  on  kz  are  chosen  as  such  since  for  |kz|>k0  the 
angles  of  propagation  are  imaginary.  Therefore,  the  source 

k 

function  is  tapered  within  the  limits  -I  < 1  . 

k0 


Finally,  to  account  for  the  half-mesh  symmetry  of  the 


..  Az 

ik,TT 


depth  grid  a  phase  term  in  the  wavenumber  domain  of  e  z 2  is 


added.  The  wavenumber  domain  starting  field  for  the  wide- 
angle  point  source  can  then  be  written  as 


T(r  =  OX)  =  -2i 


jRo 

'2nk, 


sin(k  z  ) 


l-—. 

ko; 


k2V1/4 


(3.39) 


Following  Tappert's  procedure,  for  the  COIPE  it  is  eas¬ 
iest  to  generate  the  starting  field  in  the  tilde  domain.  The 
starting  depth,  zs,  must  then  be  transformed  to  tilde  space. 

Defining 


f  =  Zs~Zj~1 
zj"zj-i 


(3.40) 


the  source  depth  in  the  tilde-domain  is  then 


zs  =  Zj _  i  +  |)f  p  Zjy= 


(3.41) 


where 

Uj  =  fCj  +  ( 1  —  f)Cj _ ,  .  (3.42) 

Note  that  Eq.  (3.42)  performs  a  linear  interpolation  of  c  at 
the  source  depth.  Since,  from  Eq.  (3.22) 
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Zj  Zj-i  -  2jr0-=± 


Zj  Zj  _  1 


■Ft +J°Ti' 


(3.43) 


then  Eq .  (3.41)  becomes 


:  ^(zj-zj-Of  - 

zs  “  zj-i+  F  ) -  “  zJ-i  + 


7uj+7^ 


(3.44) 


or 


(3.45) 


which  is  merely  Eq.  (3.19)  rewritten  using  cs  =  c(zs)  .  Also 
from  the  expression  for  the  tilde  domain  transformation. 


it  follows  that 


Therefore, 


so 


Z 

z  =  J*Vn(z',  r)dz' , 

o 


*  =  (n(z,r))"1/2 
dz  ’ 


(3.46) 


(3.47) 


(3.48) 


z  =  J(n(z,  r))-1/2dz’ .  (3.49) 

o 
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Since  c(z)  =  c(z)  ,  then  n(z)  =  n(z) ,  such  that 


2  Z 

z=  |*(n(z', r))~1/2dz’  =  -j==[ .  (3.50) 

J  Jc0(r)J 

0  0 

Equation  (3.50)  could  be  solved  numerically,  but  it  does  not 
have  the  analytical  expression  the  previous  transformation 
did.  Rather,  the  previous  transformation  for  source  depth 
can  be  applied  to  receiver  depths  at  the  original  gridding  of 
the  environment.  In  other  words,  consider 

Zk  =  (k-C)AZ  (3.51) 

which  is  essentially  what  was  obtained  in  Eq.  (3.19),  where 
now  the  tilde  depths  are  those  found  in  the  original  trans¬ 
formation  vector,  Zj  for  j  =  2...Ns-l  .  Below  Zns-i  /  Zk>Ns-i  is 
used. 

To  obtain  the  field,  however,  it  must  be  noted  that 

^(z,  r)  =  [n(z,  r)]-1/4lP(z,  r)  (3.52) 

or,  since  n(z)  =  n(z)  , 

[n(z,r)]4Y(z,r)  =  W(r)  .  (3.53) 

In  either  case,  'F(zj)  must  first  be  interpolated  from  'F(zjc)  . 

Since  there  is  no  simple  gradient  structure  for  *P(z) ,  a 
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numerical  scheme  must  be  used.  For  this  case,  a  typical  5- 


point  interpolation  scheme  is  used  to  obtain  'F(zj)  .  If  Eq. 
(3.52)  is  used,  then  n(zk,  r)  is  evaluated  and  the  values  of 

T(zk,r)  =  [n(zk,  r)]4xF(zj)  (3.54) 

are  computed  where  k  =  j-1  and  c(zk)  is  evaluated  from  the 

original  profile.  Note  that  [n(zk, r)]4'J/(Z(C)  could  have  been 
computed  and  then  interpolated  to  Zj  .  However,  this  would 

only  be  an  approximate  interpolation  of  n  .  Therefore,  the 
former  method  is  an  exact  interpolation  of  n  and  only 

approximate  in  'F ,  which  is  expected  to  be  the  more  accurate 
method . 

As  was  briefly  discussed  in  the  previous  chapter,  the 
COIPE  approach  is  expected  to  remove  the  reference  sound 
speed  ambiguity  present  in  the  WAPE.  Removing  this  sensi¬ 
tivity  should  provide  more  accurate  results  without  the  need 
to  determine  an  appropriate  reference  sound  speed  for  the 
specific  problem.  In  the  following  chapter  the  COIPE  is 
benchmarked  against  the  SPE,  WAPE,  and  normal  mode  theory. 
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IV.  NUMERICAL  RESULTS 


The  previous  chapters  present  an  outline  of  the  histor¬ 
ical  and  theoretical  background  of  the  parabolic  equation 
approximation  to  the  acoustic  wave  equation  as  well  as  the 
COIPE  implementation.  In  this  chapter,  the  focus  shifts  to 
the  analysis  of  attempts  to  improve  the  accuracy  and  valid¬ 
ity  of  the  WAPE  approximation.  First,  the  effects  of  improv¬ 
ing  the  numerical  accuracy  of  the  WAPE,  as  implemented  by  the 
MMPE  model,  is  evaluated.  Second,  the  COIPE  implementation 
is  evaluated  under  a  number  of  test  cases.  Finally,  the  SPE, 
MMPE,  and  COIPE  are  subjected  to  the  normal  mode  decomposi¬ 
tion  technique  previously  discussed. 

A.  DOUBLE-PRECISION  IMPLEMENTATION  OF  THE  MMPE 

While  the  MMPE  has  been  used  successfully  under  a  wide 
variety  of  test  cases  and  environments  for  a  number  of  years, 
it  was  recently  tested  extensively  on  a  number  of  specific 
problems.  These  tests  were  conducted  as  part  of  the  Shallow 
Water  Acoustic  Modeling  (SWAM)  Workshop  held  at  the  Naval 
Postgraduate  School  in  1999 .  The  SWAM  workshop  provided  a 
set  of  detailed  environmental  cases  that  served  as  the  basis 
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for  comparison  between  the  numerous  modeling  efforts  being 
undertaken  in  the  United  States  and  internationally.  As 
part  of  this  workshop,  the  convergence  sensitivity  of  the 
MMPE  to  the  choice  of  the  range  step,  Ar ,  was  analyzed 
(Smith,  2000) .  As  a  result  of  this  analysis,  which  is  dis¬ 
cussed  in  more  detail  later  in  this  section,  an  effort  to 
develop  a  more  numerically  precise  MMPE  implementation  was 
undertaken.  Increasing  the  numerical  precision  of  the  model 
would  also  help  answer  whether  truncation  of  the  field,  as 
postulated  in  the  introduction,  is  a  source  of  error  in  the 
model . 

In  an  effort  to  avoid  introducing  new  errors  or  anoma¬ 
lies  into  the  MMPE,  no  algorithm  changes  were  made  to  the 
program  code  itself.  Rather,  the  variables  and  array  decla¬ 
rations,  as  well  as  the  corresponding  post  processing  rou¬ 
tines,  were  recast  as  double-precision.  The  original  MMPE 
code  was,  with  the  exception  of  the  FFT  subroutine,  predomi¬ 
nantly  single-precision.  One  of  the  concerns  raised  by  this 
procedure  was  that  computational  times  would  increase  dra¬ 
matically.  This,  however,  was  not  the  case.  Since  the  bulk 
of  the  computational  work  is  done  by  the  FFT  subroutine, 
which  was  already  in  double-precision  format,  computational 
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times  only  increased  by  three  to  five  percent.  These 
increases  are  deemed  relatively  insignificant  and  do  not  in 
and  of  themselves  hinder  the  use  of  the  double-precision 
MMPE  (DP-MMPE) . 

As  a  method  of  comparison  with  the  SPE  and  existing 
MMPE,  the  DP-MMPE  was  tested  using  one  of  the  SWAM  test 
cases.  The  test  case  chosen  was  referred  to  as  SWAM  Flat-a. 
The  environment  has  an  isospeed  water  column,  c0  =  1500  m/s 

p0  =  1.00  g/cm3,  overlying  a  flat  bottom.  The  bottom  proper¬ 
ties  are  defined  every  2  km  and  are  linearly  interpolated  out 
to  a  range  of  20  km.  Compressional  attenuation  is  fixed  at 
0.1  dB/X  and  there  is  no  shear  in  the  bottom.  The  source  is 
at  a  depth  of  30  m  and  a  frequency  of  250  Hz  CW.  The  envi¬ 
ronment  is  shown  below  in  Fig.  4.1.  Fig.  4.2  shows  the  full, 
CW  field  at  250  Hz  for  the  source  at  30  m  as  obtained  from 
using  the  MMPE  model  solution. 

During  the  theoretical  development  of  the  wide-angle  PE 
approximation  it  was  noted  that  the  range  step  size  Ar ,  as 
based  on  the  centered-step  scheme,  is  roughly  third  order 
accurate.  The  implications  of  this  fact  are  that  as  the 
range  step  size  is  decreased,  the  solution  accuracy  should 
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improve  until  a  stable  solution  converges.  In  fact,  most  PE 
models  do  follow  this  behavior.  The  PE/SSF  does  not. 


ftnge(kjTf 

Figure  4.1.  The  Flat-a  test  case  environment . 

To  illustrate  the  relation  between  the  range-step  size 
and  the  covergence  of  the  solution,  a  number  of  solutions 
were  generated  for  the  Flat-a  case  discussed  above.  The  only- 
parameter  varied  in  any  of  the  cases  was  the  size  of  the 
range-step.  In  each  of  the  cases  the  transmission  loss  (TL) 
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Figure  4.2.  MMPE  field  solution  for  Flat-a  case. 

for  a  receiver  at  35  m  is  shown  in  Fig.  4.3.  The  upper  and 
lower  panels  of  Figs.  4.3  and  4.4  show  the  first  and  last  5 
km  of  the  solution,  respectively.  Note  that  the  solution 
appears  to  be  converging  as  Ar-5  m.  As  can  be  seen  in  Fig. 
4.4,  however,  the  solution  diverges  for  further  decrease  of 
Ar ;  the  solution  degrades  for  values  of  Ar  less  than  about 
5  m.  The  solution  appears  to  reach  convergence  at  Ar-A,. 
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A  possible  cause  of  this  divergence  was  postulated  to 
be  truncation  of  the  field  and  a  lack  of  numerical  precision 
in  the  declaration  of  variables  in  the  MMPE  computer  code . 
In  an  effort  to  test  this  theory  the  DP-MMPE  model  was  run 
with  identical  model  inputs  to  those  used  in  the  convergence 
and  stability  testing  shown  above  for  the  MMPE.  The  results 
of  this  testing  are  shown  below  in  Figs.  4.5  and  4.6. 


Figure  4.5.  DP-MMPE  convergence  testing  for  various  range- 

step  sizes. 


Figure  4.6.  DP-MMPE  stability  testing  for  various  range-step 

sizes . 

Clearly,  there  does  not  appear  to  be  any  discernible 
improvement  in  the  convergence  of  the  solution  through  the 
use  of  the  DP-MMPE.  It  exhibits  the  same  general  behavior  as 
that  of  the  original  MMPE  model.  A  comparison  between  the 
MMPE  and  DP-MMPE  at  the  5  m  and  2  m  range-step  sizes  is  shown 
below  in  Fig.  4.7. 
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Figure  4.7.  DP-MMPE  and  MMPE  comparison  for  range-step  sizes 

of  5  m  and  2  m. 

Referring  to  Fig.  4.7,  it  can  be  seen  that  for  a  range 
step  size  of  5  m  the  solution  produced  by  the  two  models  are 
indeed  overlapping.  For  a  range-step  size  of  2  m  the  two 
solutions  differ  slightly  with  range,  with  the  DP-MMPE 
appearing  to  diverge  less  at  some  range  points.  One  may  con¬ 
clude,  however,  that  while  the  DP-MMPE  does  no  worse,  it  does 
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not  offer  a  solution  as  to  why  the  PE/SSF  algorithm  seems  to 
have  a  lower  limit  on  the  size  of  the  range-step. 

As  a  result  of  separate  analysis  (Smith,  2000) ,  the 
primary  cause  of  this  non- convergence  was  found  to  be 
related  to  the  structure  of  the  propagator  functions . 

Recall  that  the  environmental  propagator  function  is  of  the 

form  e-lkoArUoP(r.  z)  /  and  the  wavenumber  propagator  function  is  of 

the  form  e~lkoArTop(kz)  .  The  wavenumber  propagator  decays  expo¬ 
nentially  beyond  kz>k0.  From  this  analysis  it  can  be  shown 

that  the  PE/SSF  algorithm,  for  large  range-step  size, 
attempts  to  include  a  large  amount  of  phase  information  into 
each  range-step.  If  Ar  is  chosen  to  be  too  large,  then 
errors  are  generated  in  the  solution.  Conversely,  as  Ar  is 
reduced  in  size  there  is  inadequate  phase  information  con¬ 
tained  in  each  range- step  and  the  solution  once  again 
becomes  inaccurate.  The  most  accurate  solutions  may  be 
expected  for  range-steps  where  a  full  cycle  of  phase  infor¬ 
mation  is  contained  in  each  propagator.  The  level  at  which 
Ar  appears  to  provide  the  greatest  degree  of  convergence  is 
Ar-A,  (Smith  2000)  for  shallow  water  environments  such  as 
tested  here. 
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Before  leaving  the  area  of  solution  convergence  based 
on  the  range-step  size,  a  final  comparison  between  the  MMPE 
and  the  COIPE  is  made.  While  the  COIPE  model  was  not  imple¬ 
mented  specifically  to  address  this  problem,  it  is  of  inter¬ 
est  as  to  whether  any  of  the  changes  made  in  the  COIPE  affect 
the  convergence  due  to  Ar .  Below,  in  Fig.  4.8,  the  COIPE 
model  is  compared  against  the  MMPE  for  the  same  test  case 
used  above.  The  MMPE  results  were  chosen  over  the  DP-MMPE 
since  the  DP-MMPE  solutions  did  not  offer  any  definitive 
advantages  over  the  normal  implementation.  As  Fig.  4.8 
shows,  the  COIPE  does  not  appear  to  offer  any  greater  degree 
of  convergence  for  the  given  Ar  sizes  in  this  environment. 
These  results  are  inconclusive  in  regards  to  the  COIPE 
model,  but  the  implementation  was  not  expected  to  affect  the 
Ar  sensitivity. 

Finally,  as  a  means  of  comparison  against  another  mod¬ 
eling  technique  for  this  test  case,  the  MMPE  results  for 
range- step  size  5  m  are  compared  to  those  developed  by  Mikhin 
(2000)  for  SWAM.  These  results  were  chosen  as  a  benchmark 
due  to  the  high  degree  of  agreement  with  several  test  cases 
and  several  modeling  techniques . 
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Figure  4.8.  COIPE  and  MMPE  convergence  testing. 

The  comparison  of  the  MMPE  solution  with  the  solution 
provided  by  Mikhin  is  shown  below  in  Fig.  4.9.  Solutions  for 
rang-step  sizes  of  20,  5,  and  2  m  were  chosen  to  demonstrate 
the  degree  of  convergence  with  the  reference  solution  for 
varying  Ar  .  The  figure  shows  very  good  agreement  between 
the  MMPE  and  Mikhin  solutions  at  short  ranges .  The  two  tech- 
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niques  show  a  greater  degree  of  separation  with  increasing 
range.  This  is  most  likely  the  result  of  cumulative  phase 
errors  inherent  to  the  wide-angle  PE  technique  employed 
here.  The  figure  also  demonstrates,  once  again,  that  the 
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Figure  4.9.  MMPE  and  Mikhin  solution  comparison  for  various 

range  step-sizes. 

choice  of  a  5  m  range-step  size  provide  the  greatest  degree 
of  agreement  with  the  reference  solution. 


B.  THE  COIPE  IMPLEMENTATION  OF  THE  MMPE 


As  was  discussed  in  the  last  chapter,  the  COIPE  model 
was  developed  by  Tappert  (1995)  in  an  attempt  to  remove  or 
reduce  the  dependence  of  the  wide-angle  PE  approximation  on 
the  choice  of  the  reference  sound  speed  value  while  increas¬ 
ing  the  overall  accuracy  of  the  model.  While  this  sensitiv¬ 
ity  has  been  observed  in  a  number  of  test  cases,  it  was 
extensively  studied  during  one  of  the  PE  Workshop  II  (PE-II) 
test  cases  (Porter  and  Jensen,  1993)  .  This  case,  referred  to 
as  the  Porter  duct  problem,  was  based  on  the  observation  that 
many  of  the  wide-angle  split-step  PE  approximations  overes¬ 
timated  the  value  of  the  TL  when  applied  to  long-range  prop¬ 
agation  in  a  leaky  surface  duct  environment .  The  SPE 
implementation  of  the  PE/SSF  algorithm  does  not  generate 
this  anomaly.  Additionally,  other  PE  implementations  based 
on  Pade'  series  expansions  of  the  square  root  operator 
(which  maintain  the  proper  normal  mode  basis  set)  do  not 
result  in  this  overestimation  of  the  downrange  TL.  The 
results  of  PE-II  showed  that  many  wide-angle  models  that  do 
not  use  the  Thoms on -Chapman  (TC)  wide-angle  approximation 
avoid  this  problem. 
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The  Porter  duct  environment  is  one  of  a  source  and 


receiver  located  in  a  250  m  deep  surface  duct.  The  source  is 
at  a  depth  of  25  m  and  frequency  of  80  Hz.  The  receiver  is 
at  a  depth  of  100  m.  The  ocean  bottom  is  range- independent 
with  a  lossy  fluid  half-space  beginning  at  a  depth  of  4  km. 

Of  the  78  propagating  modes  in  the  water  column,  one  is 
trapped  in  the  surface  duct.  Figures  4.10  and  4.11  illus¬ 
trate  the  sound  speed  profile  and  TL  field  for  this  environ¬ 
ment  respectively. 


Figure  4.10.  Porter  duct  sound  speed  profile  plot. 
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Figure  4.11.  TL  field  plot  for  the  Porter  duct  problem. 

Finn  Jensen,  during  PE  Workshop  II,  showed  that  the  SPE 
was  in  excellent  agreement  with  the  normal  mode  reference 
solution  used  in  the  test  case  (Porter  and  Jensen,  1993) . 
Because  this  reference  solution  is  no  longer  available,  the 
SPE  results  were  used  as  the  reference  for  this  analysis. 
The  only  alterations  to  the  various  model  input  parameters 
changed  in  the  following  analysis  is  that  of  the  reference 
sound  speed,  c0 .  Prior  to  looking  at  the  COIPE  model  or  the 
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effects  of  varying  c0,  the  TL  plot  for  the  receiver  depth  of 
100  m  is  shown  in  Fig.  4.12  for  the  SPE  and  WAPE  implementa¬ 
tions  when  a  reference  sound  speed  of  c0  =  1500  m/s  is  used. 

It  can  be  clearly  seen  that  the  WAPE  solution  shows  a  dis¬ 
tinct  drop  off  at  ranges  beyond  about  60  km.  It  is  this 
anomaly  that  the  COIPE  model  is  meant  to  eliminate. 


Figure  4.12.  TL  Plots  for  SPE  and  WAPE  at  the  receiver  depth 

of  100  m  and  c0  =  1500  m/s. 

The  reference  sound  speed  was  then  set  at  1485  m/s  and 
the  two  models  were  run  again.  These  results  are  shown  in 
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Fig.  4.13.  This  plot  illustrates  the  effect  the  value  of  c0 
has  on  the  WAPE  solution.  By  "correctly"  choosing  the  value 
of  c0,  the  WAPE  correctly  determines  the  long  range  TL.  The 
difficulty,  however,  is  in  correctly  determining  this  value. 
A  change  in  c0  by  as  little  as  5  m/s 

to  either  side  of  this  "correct"  value  results  in  the  long 
range  drop  off  shown  in  Fig.  4.12. 


Figure  4.13.  TL  Plots  for  SPE  and  WAPE  at  the  receiver  depth 

of  100  m  and  c0  =  1485  m/s. 
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The  COIPE  was  then  tested  using  identical  environmental 


inputs  to  those  used  in  the  analysis  above.  COIPE  model  runs 
were  made  at  c0  values  of  1500  and  1485  to  provide  comparison 
with  the  SPE  results.  These  are  shown  in  Figs.  4.14  and  4.15 
respectively.  The  COIPE  model  compares  well  with  the  SPE 
solution  at  all  ranges  and  there  is  no  drop  off  at  extended 
range  as  was  observed  with  the  WAPE.  The  COIPE  solution  is 
indeed  insensitive  to  changes  in  the  reference  sound 
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Figure  4.14.  TL  Plots  for  SPE  and  COIPE  at  the  receiver  depth 

of  100  m  and  cn  =  1500  m/s. 


Figure  4.15.  TL  Plots  for  SPE  and  COIPE  at  the  receiver 
depth  of  100  m  and  c0  =  1485  m/s. 

speed,  as  shown  in  Fig.  4.15.  While  both  plots  in  Fig.  4.15 
do  show  phase  differences  from  those  of  Fig.  4.14,  the  COIPE 
follows  the  reference  solution  equally  well  in  both  cases. 
The  COIPE  was  tested  at  a  wide  variety  of  c0  values  with  sim¬ 
ilarly  successful  results. 

Clearly  the  COIPE  corrects  the  TL  drop-off  encountered 
in  the  MMPE  and  similar  wide-angle  implementations.  The 


68 


COIPE  corrects  the  small  error  in  the  phase  calculation 
caused  by  the  incorrect  choice  of  the  reference  sound  speed. 

A  small  error  in  this  phase  calculation  can  cause  a  large,  on 
the  order  of  20  dB,  downrange  TL  error  (Porter  and  Jensen, 
1993)  .  The  test  case  highlights  this  problem  because  it  sup¬ 
ports  a  single  trapped  mode  in  the  surface  duct  that  is  a 
leaky  mode.  The  trapped  mode  continually  loses  energy  into 
the  lower  region  where,  due  to  the  strong  upward  refracting 
profile,  it  is  directed  back  up  toward  the  surface  duct  caus¬ 
ing  interferences.  Small  phase  errors  result  in  incorrectly 
applying  these  interferences  to  the  surface  duct  mode.  The 
transformation  and  subsequent  calculation  of  the  reference 
sound  speed  at  each  range  step  removes  the  source  of  these 
phase  errors  and  thus  correctly  applies  the  interferences 
described  above. 

C.  NORMAL  MODE  DECOMPOSITION  OF  THE  PE  FIELD 

The  mode  functions  of  the  wide-angle  PE  approximation 
form  a  different  basis  set  for  modal  expansion  than  those 
obtained  using  standard  normal  mode  theory.  As  discussed  in 
the  theoretical  development  presented  earlier,  the  SPE  does 
reproduce  the  standard  basis  set.  This  section  presents  an 
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analysis  and  comparison  of  the  modal  decomposition  of  the 
SPE,  WAPE,  and  COIPE  models. 

The  ocean  environment  chosen  for  this  analysis  is  the 
Munk  canonical  profile  (Munk,  1974) .  This  is  a  deep  ocean, 
range -independent  environment  with  sound  speed  profile  char¬ 
acterized  by 

[c(z)-c0]/c0  =  eCen-n-l),  (4.1) 

where  e  =  0.0057,  the  scaled  depth  variable  r\  =  2(z-zaxis)/B  , 

and  the  reference  sound  speed  is  1490  m/s.  The  axis  depth 
was  set  at  zaxis  =  1000m,  and  the  bottom  depth  at  4500  m.  The 

sound  speed  profile  is  shown  below  in  Fig.  4.16.  The  source 
was  placed  at  a  depth  of  1000  m  with  a  frequency  of  100  Hz. 

A  computational  range  of  100  km  was  used  for  PE  field  calcu¬ 
lations.  A  bottom  depth  of  4500  m  and  a  computational  depth 
of  8000  m  were  used. 

The  acoustic  pressure  field  was  calculated  using  each 
of  the  three  models.  The  field  is  shown  in  Fig.  4.17.  The 
resulting  fields  were  then  decomposed  into  normal  modes  and 
the  mode  amplitudes  determined  using  the  techniques  outlined 
in  Chapter  II  for  each  of  the  models.  Since  the 
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Figure  4.16.  Munk  canonical  sound  speed  profile, 
environment  chosen  is  range- independent ,  the  mode  amplitudes 
should  remain  constant  with  range.  The  amplitudes  for 
selected  modes  are  shown  in  Figs.  4.18  thru  4.20  for  the  SPE, 
WAPE,  and  COIPE  respectively.  Modes  were  chosen  to  avoid 
overlapping  plots  to  improve  readability  and  to  avoid  bottom 
interactions.  The  reference  sound  speed  is  1500  m/s  for  this 
series  of  plots.  Each  of  the  three  models  produce  satisfac¬ 
tory  results  for  the  lower  modes.  Note  the  nearly  constant 
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Figure  4.17.  TL  Field  plot  for  Munk  canonical  profile, 
amplitude  with  range  for  all  three  models.  The  higher  order 
modes  of  the  WAPE  plot,  however,  show  increasing  levels  of 
fluctuations.  Conversely,  the  COIPE  modal  amplitudes  com¬ 
pare  well  with  the  SPE  results  and  do  not  show  the  high 
degree  of  fluctuation  seen  in  the  WAPE  plots. 

The  three  models  were  then  tested  for  this  environment 
using  a  reference  sound  speed  of  1490  m/s.  The  results  are 
shown  in  Figs.  4.21  through  4.23.  The  results  for  the  COIPE, 
Fig.  4.23,  are  visibly  unchanged,  as  expected.  The  WAPE 
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amplitudes  show  somewhat  less  fluctuation  for  the  higher 
modes.  This  can  most  likely  be  attributed  to  the  sensitivity 
to  the  reference  sound  speed  discussed  in  the  previous  sec¬ 
tion. 

This  analysis  has  demonstrated  the  wide-angle  PE's  lack 
of  the  ability  to  reproduce  the  standard  normal  mode  basis 
set  as  defined  by  normal  mode  theory  and  as  duplicated  by  the 
SPE.  The  WAPE  shows  increasing  levels  of  fluctuation  in  the 
modal  amplitudes  with  increasing  mode  numbers.  This  fluctu¬ 
ation  results  in  the  produced  PE  field  not  exactly  repre¬ 
senting  the  environment  under  study.  The  COIPE,  however, 
does  not  show  the  same  degree  of  fluctuation  in  modal  ampli¬ 
tudes,  and  therefore  should  better  represent  the  acoustic 
field.  It  should  be  noted  that  the  modal  amplitudes  of  the 
COIPE  do  not  exactly  match  those  of  the  SPE.  This  is  due  to 
the  generation  of  the  source  function  in  the  tilde-domain.  A 
more  proper  treatment  might  generate  the  starting  field  in 
real  space  and  then  transform  to  tilde  space.  However,  this 
would  introduce  some  numerical  errors  to  the  transformation 
of  the  filed  form  tilde  to  real  space.  This  is  an  area  for 
further  analysis  and  study. 
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V.  CONCLUSIONS  AND  SUMMARY 


This  thesis  has  examined  the  parabolic  equation  approx¬ 
imation  to  the  acoustic  wave  equation  and,  in  particular, 
two  possible  implementations  aimed  at  removing  inherent  WAPE 
errors  were  analyzed.  The  theoretical  background  underlying 
each  of  these  models  was  presented  along  with  the  specific 
numerical  implementation  used  for  testing  and  comparison 
with  existing  PE  approximations.  Finally,  numerical  results 
were  provided  comparing  the  SPE,  WAPE,  DP-WAPE,  and  COIPE 
models  against  a  number  of  test  cases.  These  included 
results  from  the  SWAM  workshop,  the  Porter  duct  problem,  and 
the  Munk  canonical  profile. 

The  parabolic  equation  and  many  of  its  derivatives  suf¬ 
fer  from  a  number  of  anomalies  or  errors  depending  on  the 
particular  implementation  being  used.  Of  these,  three  were 
examined  in  this  analysis.  First,  the  selection  of  the 
range-step  size  used  by  the  PE/SSF  algorithm  affects  both 
the  convergence  and  the  stability  (Smith,  2000)  of  the  solu¬ 
tion  with  an  error  to  second  or  third  order  in  the  range-step 
increment  (Jensen,  et.  al.,  1994).  Second,  the  WAPE  suffers 
an  anomaly  that  under  certain  environmental  conditions  it 
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incorrectly  computes  the  down- range  TL  (Porter  and  Jensen, 
1993) .  This  anomaly  is  manifest  in  a  distinct  drop  in  TL. 
Results  of  the  Parabolic  Equation  Workshop  II  (PE- II) 
attribute  this  error  to  incorrectly  handling  the  loss  of 
energy  from  a  leaky  surface  duct  mode  and  its  subsequent 
interaction  with  the  remaining  sound  field.  It  was  shown 
here,  and  in  PE-II,  that  the  manifestation  of  this  anomaly  is 
highly  dependent  on  the  choice  of  the  reference  sound  speed. 

A  small  change  in  the  value  of  c0  results  in  a  significant 
change  in  the  computed  down-range  TL.  Third,  the  WAPE  does 
not  reproduce  the  standard  normal  mode  basis  set  as  defined 
by  normal  mode  theory  and  as  reproduced  by  the  SPE.  In  addi¬ 
tion  to  the  phase  errors  experienced  to  some  degree  by  all  PE 
approximations,  the  WAPE  mode  amplitudes  exhibit  a  range 
dependence  in  a  range- independent  environment  (Smith  and 
Smith,  1997) .  This  effect  was  found  to  contribute  to  the 
anomalous  TL  findings  just  described. 

A  double-precision  version  of  the  WAPE  model  was  imple¬ 
mented  in  an  attempt  to  remove  or  improve  the  convergence 
sensitivity  to  the  selection  of  the  range-step  size  used  by 
the  PE/SSF  algorithm.  The  double-precision  model  was  tested 
against  the  existing  SPE,  WAPE,  and  COIPE  models  in  a  simple 
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range-dependent  environment  originally  defined  for  SWAM. 

The  double-precision  model  was  found  to  do  as  well  as,  but  no 
better  than,  the  existing  approximations  under  the  tested 
conditions.  From  this  analysis  it  was  concluded  that  the 
range-step  convergence  sensitivity  is  inherent  to  the  PE/SSF 
algorithm  used  by  the  models.  This  is  supported  by  separate 
analysis  done  by  Smith  (2000) .  The  most  accurate  solutions 
may  be  expected  for  range-steps  where  a  full  cycle  of  phase 
information  is  contained  in  each  of  the  propagators  used  by 
the  algorithm.  The  level  at  which  Ar  appears  to  provide  the 
greatest  degree  of  convergence  is  Ar~A  (Smith  2000)  for 
shallow  water  environments  such  as  those  tested.  The  COIPE 
model  was  also  tested  for  sensitivity  to  the  choice  of  the 
range-step  on  solution  convergence.  It  was  found  to  be  sim¬ 
ilarly  affected,  as  was  expected. 

The  COIPE  model,  a  WAPE  implementation,  was  developed 

to  remove  or  reduce  the  solution  sensitivity  to  the  choice  of 

the  reference  sound  speed  value  ( Tapper t,  1991) .  The  COIPE 

model  uses  an  operator  transformation  which  results  in  a 

range-dependent  series  of  values  for  the  reference  sound 

speed.  The  COIPE  model  was  compared  with  the  SPE  and  WAPE 

models  using  the  Porter  duct  problem  developed  for  PE-II 
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(Porter  and  Jensen,  1993)  .  These  tests  were  conducted  using 
a  variety  of  reference  sound  speed  values,  with  successful 
results  in  all  cases.  The  COIPE  model  does  not  exhibit  the 
TL  drop  with  range  seen  with  the  WAPE. 

Finally,  the  COIPE  was  tested  against  the  SPE  and  WAPE 
models  for  its  ability  to  decompose  properly  into  the  stan¬ 
dard  mode  basis  set.  The  models  were  tested  in  a  deep-ocean 
Munk  canonical  profile  environment.  As  discussed,  the  WAPE 
does  not  correctly  decompose  with  this  basis  set,  and  shows  a 
range  dependence  in  its  modal  amplitudes  in  a  range- indepen¬ 
dent  environment.  This  range  dependence  can  be  seen  as  fluc¬ 
tuations  in  the  modal  amplitudes  with  range.  The  COIPE  modal 
amplitudes  show  much  less  fluctuation  with  range  than  those 
of  the  WAPE,  and  more  closely  reproduce  those  obtained  by  the 
SPE.  This  is  a  result  of  the  improved  handling  of  the  phase 
by  the  COIPE  as  outlined  above  and  demonstrated  through  the 
modal  decomposition. 

The  double-precision  model  testing  provides  further 
evidence,  in  addition  to  that  shown  by  Smith  (2000)  ,  that  the 
small  range-step  sensitivity  is  a  fundamental  property  of 
the  PE/SSF  algorithm,  and  no  further  investigation  of  this 
property  is  necessary.  The  results  of  the  COIPE  model  test- 
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ing,  while  promising,  are  not  conclusive.  The  COIPE  cor¬ 
rects  the  error  seen  in  the  Porter  duct  problem,  and  improves 
the  reproduction  of  the  standard  normal  mode  basis  set,  but 
requires  further  investigation.  Preliminary  testing  in 
shallow  water  environments  indicate  the  COIPE  may  not  pro¬ 
duce  results  significantly  different  from  those  of  the  WAPE 
model.  Additionally,  the  mode  amplitudes  themselves,  while 
showing  less  fluctuations,  are  of  different  relative  magni¬ 
tudes  from  those  produced  by  the  SPE.  This  remains  as  an 
area  for  further  investigation. 


81 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


82 


LIST  OF  REFERENCES 

Bellman,  R. ,  (1964) .  Perturbation  Techniques  in  Mathemat¬ 

ics,  Physics,  and  Engineering  (Holt,  Rinehart,  and  Winston, 
Inc.,  New  York),  pp.  38-39. 

Collins,  M.  D.,  (1994).  "Generalization  of  the  Split-Step 

Pade  Solution,"  J.  Acoust.  Soc.  Am.  96,  pp.  382-385. 

Fock,  V.  A. ,  (1965) .  Electromagnetic  Diffraction  and  Propa¬ 

gation  Problems  (Pergamon,  New  York) . 

Hardin,  R.  H.  and  Tappert,  F.  D.,  (1973).  "Applications  of 

the  Split-Step  Fourier  Method  to  the  Numerical  Solution  of 
Nonlinear  and  Variable  Coefficient  Wave  Equations,"  SIAM 
Rev.  15. 

Jensen,  F.  B.,  Kuperman,  W.  A.,  Porter,  M.  B. ,  and  Schmidt, 
H.,  (1994).  Computational  Ocean  Acoustics  (AIP,  New  York). 

Lee,  D.,  Botseas,  G.  and  Papadakis,  J.  S.,  (1981).  "Finite- 

Difference  Solution  to  the  Parabolic  Wave  Equation,"  J. 
Acoust.  Soc.  Am.  70,  pp .  795-800. 

Leontovich  and  Fock,  V.  A.,  (1946).  Zh.  Eksp.  Teor .  Fiz. 

(Sov.  Phys .  JETP,  Moscow),  pp.  13-24. 

McDaniel,  S.  T.,  (1975)  .  "Propagation  of  Normal  Modes  in  the 

Parabolic  Approximation,"  J.  Acoust.  Soc.  Am.  57,  pp.  307- 
311. 

McDaniel,  S.  T.,  (1975).  "Parabolic  Approximations  for 

Underwater  Sound  Propagation,"  J.  Acoust.  Soc.  Am.  58,  pp. 
1178-1185. 

Mikhin,  D. ,  (1999).  "Energy-Conserving  and  Reciprocal  Solu¬ 

tions  for  Higher-Order  Parabolic  Equations, "  unpublished 
report  from  Shallow  Water  Acoustic  Workshop. 

Munk,  W.  H. ,  (1974).  "Sound  Channel  in  an  Exponentially 

Stratified  Ocean,  with  Application  to  SOFAR, "  J.  Acoust. 

Soc.  Am.  55,  pp.  220-226. 


83 


Porter,  M.  B.,  and  Jensen,  F.  B.,  (1993).  "Anomalous  Para¬ 

bolic  Equation  Results  for  Propagation  in  Leaky  Surface 
Ducts,"  J.  Acoust.  Soc.  Am.  94,  pp.  1510-1516. 

Smith,  A.  R.,  and  Smith,  K.  B.  ,  (1997).  "Mode  Functions  for 

the  Wide-Angle  Approximation  to  the  Parabolic  Equation,"  J. 
Acoust.  Soc.  Am.  103,  pp.  814-821. 

Smith,  K.  B.,  (2000) .  "Convergence,  Stability,  and  Vari¬ 

ability  of  Shallow  Water  Acoustic  Predictions  Using  a  Split- 
Step  Fourier  Parabolic  Equation  Model,"  unpublished  report 
from  Shallow  Water  Acoustic  Workshop. 

Smith,  K.  B.,  Brown,  M.  G. ,  and  Tappert,  F.  D. ,  (1992).  "Ray 

Chaos  in  Underwater  Acoustics,"  J.  Acoust.  Soc.  Am.  91,  pp. 
1939-1949. 

Tappert,  F.  D. ,  (1977).  "The  Parabolic  Approximation 

Method, "  in  Wave  Propagation  and  Underwater  Acoustics,  Lec¬ 
ture  Notes  in  Physics,  edited  by  J.  B.  Keller  and  J.  S. 
Papadakis  (Springer-Verlag,  Berlin),  Vol.  70,  pp.  224-287. 

Tappert,  F.  D. ,  (1991).  "Density  in  PE/SSF, "  unpublished 

lecture  notes. 

Tappert,  F.  D. ,  and  Brown,  M.  G.,  (1996).  "Asymptotic  Phase 

Errors  in  Parabolic  Approximations  to  the  One-Way  Helmholtz 
Equation,"  J.  Acoust.  Soc.  Am.  99,  pp.  1405-1414. 

Tappert,  F.  D. ,  Spiesberger,  J.  L.,  and  Boden,  L. ,  (1995). 

"New  Full-Wave  Approximation  for  Ocean  Acoustic  Travel  Time 
Predictions,"  J.  Acoust.  Soc.  Am.  97,  pp.  2771-2782. 

Thomson,  D.  J. ,  (1993).  "PE-Based  Spectral  Decomposition," 

in  PE  Workshop  II:  Proceedings  of  the  Second  Parabolic  Equa¬ 
tion  Workshop,  edited  by  S.  A.  Cin-Bing,  D.  B.  King,  J.  A. 
Davis,  and  R.  B.  Evans,  Naval  Research  Laboratory  NRL/BE/ 
7181-93-001  (US  Government  Printing  Office),  pp.  296-320. 

Thomson,  D.  J. ,  and  Bohun,  C.  S.,  (1988).  ”A  Wide-Angle  Ini¬ 

tial  Field  for  the  Parabolic  Equation  Models,"  J.  Acoust. 
Soc.  Am.  83,  pg.  S118. 


84 


Thomson,  D.  J.,  and  Chapman,  N.  R. ,  (1983).  "A  Wide-Angle 

Split-Step  Algorithm  for  the  Parabolic  Equation,"  J.  Acoust. 
Soc.  Am.  74,  pp.  1849-1854. 


85 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


86 


1 


INITIAL  DISTRIBUTION  LIST 

1.  Defense  Technical  Information  Center . 2 

8725  John  J.  Kingman  Rd.,  STE  0944 

Ft.  Belvoir,  Virginia  22060-6218 

2.  Dudley  Knox  Library . 2 

Naval  Postgraduate  School 

411  Dyer  Rd. 

Monterey,  California  93943-5101 

3.  Prof.  Kevin  B.  Smith . 3 

Department  of  Physics,  CODE  PH/SK 

Naval  Postgraduate  School 
Monterey,  California  93943-5117 

4.  Prof.  James  V.  Sanders . 1 

Department  of  Physics,  CODE  PH/SJ 

Naval  Postgraduate  School 
Monterey,  California  93943-5117 

5.  Undersea  Warfare  Programs  Office,  Code  35 . 1 

833  Dyer  Rd.,  Room  304 

Naval  Postgraduate  School 
Monterey,  California  93943-5124 

6.  Kirk  A.  Weatherly,  LT  USN . 2 

607  South  Main  Street 

Farmland,  Indiana  47340 

7.  Joseph  D.  Weatherly . 2 

St.  Rd.  32  W,  PO  Box  202 

Farmland,  IN  47340 


87 


